How to overcome the numerical instability of 
the scheme of divided differences? 

Alicja Smoktunowicz, Przemyslaw Kosowski and Iwona Wrobel 

Faculty of Mathematics and Information Science, Warsaw University of 
Technology, PI. Politechniki 1, 00-661 Warsaw, Poland 

Abstract 

The scheme of divided differences is widely used in many approx- 
imation and interpolation problems. Computing the Newton coeffi- 
cients of the interpolating polynomial is the first step of the Bjorck 
and Pereyra algorithm for solving Vandermonde systems of equations 
(Cf. pQ). Very often this algorithm produces very accurate solution. 
The problem of determining the Newton coefficients is intimately re- 
lated with the problem of evaluation the Lagrange interpolating poly- 
nomial, which can be realized by many algorithms. For these reasons 
we use the uniform approach and analyze also Aitken's algorithm of 
the evaluation of an interpolating polynomial. We propose new algo- 
rithms that are always numerically stable with respect to perturbation 
in the function values and more accurate than the Aitken's algorithm 
and the scheme of divided differences, even for complex data. 

AMS subject classification: 65D05, 65G50 
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1 Introduction 

The aim of this paper is to study the numerical stability of algorithms of eval- 
uating the Lagrange interpolating polynomial w N (z) = Y,f=o fj Hi=o,i^j J^l 
and computing the Newton coefficients c , . . . , of w^{z). We assume that 
z , zi, . . . , z N are pairwise distinct data points and /o, /i, • • • , /jv are associ- 
ated values (real or complex) of function /. Then w n(z) interpolates function 
f(z) at points zq, z±, . . . , zn, i-e. f(zj) = fj = WN^Zj) for j = 0, 1, . . . , N and 
^n(z) = Ef =0 c j UlZi (z - Zi). 

The problem of determining the Newton coefficients is intimately related 
with the problem of evaluation the Lagrange interpolating polynomial, which 
can be realized by many algorithms. For these reasons we use the uniform 
approach and consider the following general problem. 

Problem I. For given t G C, Zo, z±, . . . , Zn G C (distinct) and /o, /i, ■ ■ ■ , /at 
compute 

n n 

Pn(z;t) = J2f ] II T^rf' n = 0,l,...,JV. (1) 

j=0 i=0, ijtj ^ 1 

In the special case of z = tzj for some j we obtain p n (z, t) = fj t n . 
Notice that for t = 1 we have 

n n 

Pn (z;l)=w n (z) = J2fi II n = 0,l,...,iV (2) 

i=o i=o, ^ * 

and for t = 0, z = 1 we get 

n n j 

Pn(l;0) = c n = X;/ i 11 ^37' " = 0,l,...,iV. (3) 

j=0 i=0, «t^j ^ * 

Usually the Newton coefficients are computed by the scheme of divided 
differences. However, it is known (Cf. jS]) that the reputation of this algo- 
rithm is not so good. Its numerical properties depend on the distribution and 
ordering of the interpolation points. A. Kielbasihski pp. 51-53) proved 
that if the nodes are real and monotonically ordered, i.e. z < z\ < . . . < z^ 
or zq > z\ > . . . > zn, then every Newton coefficient computed in fl is the 
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exact divided difference for slightly perturbed values fj (see Proposition 1), 
that is numerically stable in the following sense. 

Definition 1 We say that an algorithm W of computing 

n n 

p n (z;t)=J2fr II 7~r7' n = 0,l,..-,N 

j=0 i=0,i^j j 1 

is numerically stable (with respect to the data f , / 1; . . . , f N ) if the values 
p n (z;t) computed by W in floating point arithmetic satisfy 

n n 

Pn(z;t) = J2m + ^ n) )] II ^ n) \^e M k N (4) 

3=0 i=0,i^j j 1 

and kjy is a modest constant. 

It is easy to check that (4) is equivalent to 

\p n (z;t) -p n (z;t)\ ^ e M k N ^2 \fj\ II j _ V n = 0,...,N. (5) 

3=0 i=0, i+j \ Zj Z ^ 

For all set of values we have 

^max^ \p n (z; t) - p n (z; t)\ 

— < e M k N cond f (z , . . . ,z N ;z;t), (6) 



max \p n (z;t)\ 

n=0,...,N 

where 



\Z — tZj 



max > /,• If - 

cond f ( Z(h ...,z N ;z;t) = max \ Pn{z - t) \ ( 7 ) 

n=0,...,N 

is the measure of sensitivity of Problem I with respect to the data /o, /i, • • • , /at 
and will be referred to as the condition number. 

In Section 1 we review two standard methods, i.e. Aitken's interpolation 
algorithm and the scheme of divided differences. Section 3 presents a new 
efficient algorithm for computing (1). Error analysis of these algorithms is 
given in Section 4. The new algorithm is shown to be numerically stable in 
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sense (4), even for complex data, opposite to the scheme of divided differ- 
ences and Aitken's algorithm. We also show that these two algorithms are 
numerically stable if the knots are real and monotonically ordered (Cf. [I], 
pp. 51-53, [Oj). Numerical tests in MATLAB included in Section 5 confirm 
theoretical results. 

The new algorithm doesn't require any special ordering of interpolation 
points, so we can try to reorder the knots to achieve small relative error (see 
(6)). However, it seems to be a hard problem to minimize a condition number 
of Problem I (see (7)). It was observed that the Leja ordering (Cf. [2], [E|) 
may be a good choice. To reorder the points z , zi, . . . , z^ by Leja, first choose 
z satisfying \zo\ = maxj \zj\, and then determine Zj, j = 1, 2, . . . , N — 1 such 

that ni=o \ Z J ~ z k\ = m & x i>j llCo z k \. 

Although in many cases helpful, sometimes it gives results comparable 
to or even worse than monotone ordering, even for equidistant points (see 
Section 5, Table 3). 

2 Algorithms 

The simplest way to compute Newton coefficients (3) is the well known 
scheme of divided differences, which can be visualized as follows: 



\ c 
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where Cj are the divided differences of nth order. They are defined in the 
following way 

(n-l) _ (n-1) 
tt v - ~ A _ J") _ S S'+l 

J v j' i+i' • • • ' ^j+nj "-"j 

Z i ~~ z j+n 

for j = 0, 1, . . . , JV - n, and n = 1, 2, . . . , JV. 

Computing the value of the Lagrange interpolating polynomial at a given 
point z (see formula (2)) can be realized by Aitken's algorithm (see [3], 
Problem 8 on p. 289). 

The following Algorithm I can be used for computing both Newton coef- 
ficients and the value of interpolating polynomial, depending on z and t. For 
t = and z — 1 algorithm reduces to the divided differences scheme, and for 
t — 1 it is Aitken's algorithm. For all other choices of z and t Algorithm 1 is 
a general way of evaluating (1). 

Algorithm I 

For given (real or complex) z n and f n , t and z this algorithm computes 

n n 

Pn (z;t) = J2fj II J7ZJ.> n = 0,l,...,N. 

j=0 i=0, tyj 3 1 

Choose t and z 

cf :=fj, ./ 0.1 V 

for n = 1,2, ...,N 

for jfe = 0, 1, . . . , JV - n 

r„) [z - tz n+k ) 4" ^ - (z - c fe ^_ 1 1 
C * := (8) 

•°K ^k+n 

end 

p n (z;t) := c[, ra) 

end 

Then 

k+n k+n 

c i n) = E^' II n = 0,l,...,JV, k = 0,l,...,N-n. (9) 
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It means that for t = 1 c h is the Lagrange interpolating polynomial based on 
the knots z&, -Zfc+i, . . . , Zfc +ri and for z = 1, £ = is the divided difference 
of nth order based on the knots . . . , Zfc+n- 

The complexity of the divided differences and Aitken's algorithms is 
0(N 2 /2) and 0(3N 2 /2) complex multiplications, respectively. The general 
algorithm requires 0(5N 2 /2) complex multiplications. 



3 New algorithm 

In this section we propose a new algorithm for computing (1). The idea is 
similar in spirit to the ones introduced in and [Zj, but our Algorithm 

II is more efficient and robust (see Theorem 2). 
For z ^ tzj, j = 0, 1, . . . ,N 

JL a 

II ( z - tz i)= -ZT~' n = 0,l,...,N, 

2 XZj 

8=0, 17=7 J 



where 



A n = (z - £Z )(Z - ^l) • • • • • - tZ n ). (10) 



Thus 



Pw (z;o^ae , _ f \ (») » ^ n)= n ( n > 

i=o ( z w j i=o, i^j 

Note that A n has the following form 



A, 



(z - Zq) ■ . . . ■ (z - z n ) for t = 1, 

1 for t = 0, z = 1. 



We observe that 



') =1, = (%• - z n ) for j = 0, 1, . . . , n - 1, 

A n = (z- tz n )An-i, n = 0, 1, . . . , N, with A-\ = 1. 
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Then for n — 0, 1, . . . , N we have p n (z; t) = A n B n where 

n „ 

*»=£#\ = : _!\ w J = M,...,n. (12) 

j"=0 \ Z tZ 3> W j 



We see that 



and 



6 (o) = ^L_ J = 0) l,...,iV (13) 

6 (0) 

^ n) = ^y, J = 0,l,...,n. (14) 
Notice that if -B„_i is written in the form 

B n _ l = bt l) +b { r l) + ---+bti\ 

then we can rewrite B n as follows 

i/n-l) r(n-l) ,(n-l) ,(0) 

B n = -A t+t^ r+. • ■+ , ^ , 4 



(^o — Z n ) (zi — z n ) (z n _i — z n ) (z n — Zo){Z n — Zi). . \Z n — Z n -ij 

Further, the formulae (11)-(12) let us write 

b (n-i) 

6 («) = ^_ > j = 0,l,...,n-l, n = l,2,...,N (15) 

Zj z n 

and 

rij=o ( Zn ~ z i) 



W = ^ (16) 



Algorithm II 

For given (real or complex) z n and f n , t and z ^ tZj this algorithm 
computes 

n n ^ 

p n (z;t) = J2fj II J-T7' n = 0,l,...,N. 

3=0 i=0,ift j 1 

Choose t and z 



b (o) — fj — ? ' = 01 N 

Z ~CZ' 



°3 

A := z - tz 
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B := C } 

for n = 1, 2, . . . , N 
A n '.= [z — tz n )A n _i 
for j = 0, 1, . . . , n — 1 



end 



3 ? . — y 



b 



(0) 



B ■= T n b (n) 
p n {z;t) := A n B n 

end 

The complexity of Algorithm II is 0(N 2 ) flops. 

Note that this algorithm can be implemented more efficiently than it's 
written now, namely can be computed along with b^ in the inner loop 
for 

6 (»). (- 1 ) n& « ) 



llj=0 i Z 3 Z n) 

We decided to write Algorithm II this way to make it more clear. 



4 Error analysis of algorithms 

We consider complex arithmetic (cfl) (Cf. .0]) implemented using real arith- 
metic (fi) with machine unit €m- We assume that for x, y G K we have 

fl(x±y) = (x±y)(l + 6), \8\^e M . (17) 

Then 

cfl(x®y) = (xOy)(l + 5), \5\^c G e M forx,yeC, (18) 

where 

1 for 0G {+,-}, 
c = { 1 + y/2 for = *, (19) 
4 + y/2 for = /. 
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The constant c Q is not significant, it depends on the implementation of float- 
ing point operations. In fact, the bounds of errors in Algorithms I and II 
are of the same form for real and complex arithmetic, only the constants are 
increased appropriately. 

Our error analysis is uniform, it includes both real and complex cases, 
i.e. the data Zj, fj, z and t can be either real or complex. 

The values computed in fl or cfl will be marked with tilde in order to 
distinguish them from the exact ones. 

4.1 Error analysis of Algorithm I 

Theorem 1 Assume that Zj, fj,E C for j = 0,1, ... ,N and z, t G C are 
exactly representable in cfl and 

cfl(tz j )=tz ] , 3 = 0,1,...,N. (20) 

Then the values cj^ for n = 0,1, . . . N , k = 0, 1, . . . , N — n computed in cfl 
by Algorithm I satisfy 

k+n k+n 

4 n) = £^( 1 + II \Sl1\^e M ndL^ + 0(e> M ), (21) 

j=k i=k,iytj 3 



where 



, 5 in the real case, , . 

d={ (22) 
8 + 2y2 in the complex case 

and 

\Z; L Zj\ + \Zj Zk\ /oo\ 

L = max p . (23) 

0^i<j<k^N \Zi — Zk\ 

Proof. By induction on n. It follows from (8) and (17)-(19) that for n = 
0,1, ... N, k = 0,1, . . . , N — n the computed quantities satisfy 

An) _ (1 + a k n) )(z - tz n+k )ct l) - (1 + P ( k n) )(z - tz k )^ 

c k ■- _ _ _ > l 24 J 

^k ^k+n 

where 

\* k n) \,\P ( k n) \^e M d + 0(e 2 M ). (25) 
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It is clear that for n = 1 the result is immediate, since 

4" = a a + 4 1, )^ ±i + A*, (i + >)-^-, 

^fc — ^fc+1 z k+l — %k 

so we can take 5^ k = and 5^ k+1 = (3 k \ k = 0, 1, . . . , N - 1. 

Now assume that (21) is true for n — 1. We will prove that it holds also 
for n. By assumption for n — 1 we get 

fc+n— 1 fc+n 
j=fc i=fc, i^j J 



fc+n fc+n 

(* - t, fc ) = e [/.-a + oi n ^rf & - **)■ ( 27 ) 

j'=fc+l i=fc, i^j ^ 

To simplify the notation let us define 

l + ^g = (l + aj ,) )(l + ^7 1) ) 

and 

i+^S = (i+^ n) )(i+0> 

bounded as follows 

l</& l^gl < e M d (1 + (n - 1)L- 2 ) + 0(4). (28) 
Subtraction of (26) from (27) leads to 

fc+n fc+n 
(n) _ V"^ , /, r(n)s TT ^ — £2i 



where 



j=fc i=fc, i=£j Z i Z% 



1 | J(W ) _ (1 + - gfc+n) + (1 + V>g)(*fc - *j) 

%k ^fc+n 

for j = k + 1, . . . , k + n — 1, and 

1 + °fc,fc — 1 ^ 9k,ki 1 + °fc,fc+n — 1 ^ rfc,fc+ n - 

This leads to the estimation 

|5gKmax{|0g|, |^>|}L, for j = k, . . . ,k + n, 
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which, in combination with (28) and the fact that L ^ 1 (see (23)) results in 
\5^\ ^ 6MndL n ~ l + 0{ej A ), which completes the proof. ■ 

Although, in general, the constant L defined in (23) might be large, in 
the special cases the Algorithm I is numerically stable in sense (4). This hap- 
pens, for example, for specially ordered knots. For details see the following 
proposition and remarks below it. 

Proposition 1 Assume that Zj, fj G R for j — 0, 1, . . . , TV and z, t G R are 

exactly representable in fl and (20) holds. If the knots Zj are monotonically 
ordered then the values p n (z; t) = computed in fi by Algorithm I satisfy 

n n , 

p n (z;t) = J2m + ^ n) )} II \5f\^e M 5n + 0(e 2 M ). (29) 

j=0 i=0,i¥=j j 1 

Proof. It is a direct consequence of Theorem 1 for k = 0. For monotonically 
ordered knots Zj the constant L in (23) is 1. ■ 

Remark. If zj G C, zj = a + jh, j = 0, 1, . . . , N, h = (b-a)/N, a, b G C, 
then L — 1 and Algorithm I is also numerically stable, i.e. the condition (4) 
holds with the constant k N m (8 + 2y/2)N. 

In Section 5 we present numerical experiments showing that Algorithm I 
is not always numerically stable. 

4.2 Error analysis of Algorithm II 

Theorem 2 Assume that Zj, fj GC, z/ tZj for j = 0,1, ... ,N and z,t G 
C are exactly representable in cfl and (20) holds. Then the values p n (z;t) 
computed in cfl by Algorithm II satisfy 

n n 

p n (z;t) = J2lf^ + ^ n) )] II 7^5' \^ n) \^e M k N + 0(e 2 M ), (30) 
where 



3=0 i=0,i^j Zj 



5 in the real case, , , 

k N = (N + l){ (31) 
8 + 2y2 in the complex case. 
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Proof. The proof is straightforward. Using (13)-(16) and (17)-(19) we obtain 
the following quantities, computed in fl or cfl 

bf> = ^-(l + Af), j = 0,l,...,N, 

J Z — tZj J 

and for n — 1, 2, . . . , N 

A n = A n (l + On), 

&S n) = -5*— (l + Zt'), J = 0,l,...,n-1, 



Zj z n 



£(o) 



From (12) we obtain 

n 
j=0 



where 



p„(z;*) = A nJ B„(l + 5 n ), 



ia(o)| ,o(n)| < I 2 (real case) ( 2 > 

5 + V2 (complex case) 

2n + l (real case) 2 
< cm ^ +C(e M ), 
(2 + v2) n + 1 (complex case) 

In (real case) „ x 

(2 + v2) n + 3 (complex case) 

, 1 (real case) , 

1 + v2 (complex case) 



Now the formulae (15), (16) yield for all b 



(n) 
3 



If =bf\l + <pf) 1 (32) 
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where 



2 (real case) 
5 + v2 (complex case) 



+ 0(e 2 M ) 



(33) 



and for B n , n — 1,2 ... , 



iV 



5, 



j=0 



where 



3 (real case) 
6 + v^2 (complex case) 



+ 0{e 2 M ). 



Finally, this and (32), (33) lead to (30), which proves the theorem. ■ 

Remark. In the case of divided differences algorithm (formula (3)) the 
constants are smaller, namely (Cf. 0) 



5 Numerical experiments 

To illustrate our results we present numerical experiments carried out in 
MATLAB with unit roundoff e M = 2.2e-16. 

We implemented all methods and performed many tests in order to inves- 
tigate the behaviour of these methods, i.e. to compare the errors generated 
by each of them. This paragraph contains the results of some of these tests. 

Algorithm II is general, two special cases are computing the Newton co- 
efficients and computing the value of an interpolating polynomial. These are 
the cases we considered in our tests. 

The first kind of tests involves computing the Newton coefficients. One 
type of the functions we interpolated was f(z) = z s , s = 2, 3, . . . , 7. 




(real case) 
(complex case) 
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To compare our algorithm with the classical scheme of divided differences 
we evaluated the following quantities: 



err or \ 



max 

n=s+2,...,JV 



^ n ™e>£^i n 

3=0 



(34) 



and 



err or 2 = 



i=0, ij^j 



max c n 

n=s+2,...,N 

€m max IcJ 

n=0,l,...,N 



I 



(35) 



where c n and c n are the computed and the exact Newton coefficients, respec- 
tively. Notice that if f(z) = z s , then c s+2 = ■ ■ ■ = c N = for N ^ s + 2. 

For numerically stable algorithms err or 1 is small, it is a measure of insta- 
bility of an algorithm, which comes from the following chain of inequalities 
(see (4)-(7)) 



max 

n=s+2,.. 



n n 1 

|c n |^ max \c n -c n \^e M k N max V] |/_,-| If , 

,N' n=0,...,N' n=0,...,N^ J . A A . ~ 



(36) 



j=0 i=0, j^j 



Now it is easy to see that error! should be less than or equal to the constant 
k N in (4). 

Tables 1 and 2 contain the respective values of errorl and error2 for 
Aitken's algorithm and Algorithm 11 applied to the function f(z) = z 7 . The 
interpolation was based on 50 or 80 random knots ordered either monoton- 
ically or by Leja ordering. In all three tests problem was well-conditioned, 
with condition number (7) of order 1. 



Table 1: The error (34) for the scheme of divided differences (Div. diff.) and Algorithm II 
for random knots, increasing (j) and Leja orderings. 



N 


Div. diff. + t 


Alg. II + t 


Div. diff. + Leja 


Alg. II + Leja 


50 


1.161255e+04 


5.470269e-02 


4.587882e-01 


6.130860e-05 


50 


8.622460e+02 


4.401986e-02 


8.014843e-05 


3.478342e-05 


80 


2.274724e+06 


1.340639e-01 


1.195999e-04 


7.967939e-05 
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Table 2: The error (35) for the scheme of divided differences (Div. diff.) 
and Algorithm II for the same data as in Table 1. 



N 


Div. diff. + t 


Alg. II + t 


Div. diff. + Leja 


Alg. II + Leja 


50 


2.019597e+04 


9.513618e-02 


4.587882e-01 


6.130860e-05 


50 


9.928299e+02 


5.068650e-02 


8.014843e-05 


3.478342e-05 


80 


2.278485e+06 


1.342856e-01 


1.195999e-04 


7.967939e-05 



Now let us turn our attention to evaluation of an interpolating polyno- 
mial. 

Tests presented below were performed for the function f(z) = z 7 . In- 
terpolation was based on 100 equally spaced knots from the interval [—1, 1], 
ordered either monotonically or by Leja ordering. 

All graphs present the error (see (4)-(7) for the function f(z) = z 7 ) 



computed for equally spaced points from interval [—98/99,97/98], where 
p n {z, 1) is the value of an interpolating polynomial p n (z, 1) in (2) given either 
by Aitken's algorithm or by Algorithm II. 

Figures 1 and 2 describe the results for Aitken's algorithm and Algorithm 
II, respectively, for monotonically ordered knots. They confirm theoretical 
results, that Aitken's algorithm is numerically stable for monotonically or- 
dered points (increasingly or decreasingly) , see Theorem 1. 



err or 3 = 



\p n (z, 1) - z 7 



(37) 
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2.5 




Figure 1: The error (37) for Aitkcn's algorithm for f(z) = z 7 
with the knots being 100 equally spaced points 
from the interval [—98/99,97/98], increasingly ordered. 

3.5 | 1 1 1 1 1 1 1 1 1 1 




-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 



Figure 2: The error (37) for Algorithm II for f(z) = z 7 
with the knots being 100 equally spaced points 
from the interval [—98/99,97/98], increasingly ordered. 

Figures 3 and 4 present the values of error (37) for Aitken's algorithm 
and Algorithm II, respectively, for knots ordered by Leja ordering. 
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-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 



Figure 3: The error (37) for Aitken's algorithm for f(z) = z 7 
with the knots being 100 equally spaced points 
from the interval [—98/99, 97/98], ordered by Leja ordering. 



5 - 




-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 



Figure 4: The error (37) for Algorithm II for f(z) = z 7 
with the knots being 100 equally spaced points 
from the interval [—98/99, 97/98], ordered by Leja ordering. 



Table 3 contains values of error?) at five of the points in the data used to 
create figures 1-4. 
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Table 3: The error (37) for Aitken's algorithm and Algorithm II 

with knots ordered either increasingly or by Leja ordering. 



z 


Aitken + | 


Alg. II + t 


Aitken + Leja 


Alg. II + Leja 


-.989899 


8.078746e-02 


4.623302e-01 


3.791592e+06 


3.593606e-02 


-.791929 


7.973497e-02 


2.746146e-01 


1.640714e+07 


6.377801e-02 


-.593960 


3.579114e-02 


2.806026e-01 


6.421192e+05 


5.397304e-01 


-.395990 


1.215628e+00 


6.078142e-01 


7.446058e + 04 


1.823442e + 00 


-.198021 


2.787772e-01 


1.115109e+00 


3.299607e + 03 


1.115109e + 00 



Note that Leja ordering, although in many cases helpful, is not a general 
remedy. Aitken's algorithm is more accurate for monotonically ordered knots 
than for Leja points. This is because of the constant L in Theorem 1. For 
monotone order L is always 1, for Leja ordering L may be big, for example 
in the test presented above L = 197. 
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